MSE of MOMs and Relative Efficiency
Oliver de Pug
Problem 1
The method of moments estimators are not unique. Which is better?
mom1 = function(d, i){sum(d[i])/length(d[i])}
mom2 = function(d, i){n <- length(d[i])
sqrt(sum(d[i]^2)/n - mean(d[i])^2 )}
mom3 = function(d, i){sqrt(sum(d[i]^2)/(2*length(d[i])))}
x <- rexp(1000, 1/5)
mean(x)## [1] 5.381943 mom1(x)## [1] 5.381943 mom2(x)## [1] 5.31566 mom3(x)## [1] 5.348904Now bootstrap to determine the biases and standard errors.
p_load(boot)
### Use the boot function to run the bootstrap
b1 <- boot(x, mom1, R=9999)
b1##
## ORDINARY NONPARAMETRIC BOOTSTRAP
##
##
## Call:
## boot(data = x, statistic = mom1, R = 9999)
##
##
## Bootstrap Statistics :
## original bias std. error
## t1* 5.381943 0.0002329561 0.1681504 b2 <- boot(x, mom2, R=9999)
b2##
## ORDINARY NONPARAMETRIC BOOTSTRAP
##
##
## Call:
## boot(data = x, statistic = mom2, R = 9999)
##
##
## Bootstrap Statistics :
## original bias std. error
## t1* 5.31566 -0.004262346 0.2098238 b3 <- boot(x, mom3, R=9999)
b3##
## ORDINARY NONPARAMETRIC BOOTSTRAP
##
##
## Call:
## boot(data = x, statistic = mom3, R = 9999)
##
##
## Bootstrap Statistics :
## original bias std. error
## t1* 5.348904 -0.002782499 0.1749432 plot(b1)
plot(b2)
plot(b3)
x <- rexp(1000, 5)
b1 <- boot(x, mom1, R=9999)
b1##
## ORDINARY NONPARAMETRIC BOOTSTRAP
##
##
## Call:
## boot(data = x, statistic = mom1, R = 9999)
##
##
## Bootstrap Statistics :
## original bias std. error
## t1* 0.1942056 -3.612119e-05 0.006163828 b2 <- boot(x, mom2, R=9999)
b2##
## ORDINARY NONPARAMETRIC BOOTSTRAP
##
##
## Call:
## boot(data = x, statistic = mom2, R = 9999)
##
##
## Bootstrap Statistics :
## original bias std. error
## t1* 0.1938152 -0.0002113673 0.008548984 b3 <- boot(x, mom3, R=9999)
b3##
## ORDINARY NONPARAMETRIC BOOTSTRAP
##
##
## Call:
## boot(data = x, statistic = mom3, R = 9999)
##
##
## Bootstrap Statistics :
## original bias std. error
## t1* 0.1940105 -0.0001318134 0.006785682 plot(b1)
plot(b2)
plot(b3)
x <- rexp(1000, 1)
b1 <- boot(x, mom1, R=9999)
b1##
## ORDINARY NONPARAMETRIC BOOTSTRAP
##
##
## Call:
## boot(data = x, statistic = mom1, R = 9999)
##
##
## Bootstrap Statistics :
## original bias std. error
## t1* 0.9848106 0.0002115691 0.03128229 b2 <- boot(x, mom2, R=9999)
b2##
## ORDINARY NONPARAMETRIC BOOTSTRAP
##
##
## Call:
## boot(data = x, statistic = mom2, R = 9999)
##
##
## Bootstrap Statistics :
## original bias std. error
## t1* 0.984936 -0.0007004861 0.04143164 b3 <- boot(x, mom3, R=9999)
b3##
## ORDINARY NONPARAMETRIC BOOTSTRAP
##
##
## Call:
## boot(data = x, statistic = mom3, R = 9999)
##
##
## Bootstrap Statistics :
## original bias std. error
## t1* 0.9848733 -0.0003307651 0.0334574 plot(b1)
plot(b2)
plot(b3)
Problem 3
This problem deals with the properties of the MSEs. We look at their values for small and large n.
bs <- function(n, sigma2=1){0}
bsighat <- function(n, sigma2=1){sigma2/n}
vars <- function(n, sigma2=1){2*sigma2^2/(n-1)}
varsighat <- function(n, sigma2=1){2*sigma2^2*(n-1)/(n^2)}
mses <- function(n, sigma2=1){bs(n, sigma2)^2 + vars(n, sigma2)}
msesighat <- function(n, sigma2=1){bsighat(n, sigma2)^2 + varsighat(n, sigma2)}
re <- function(n, sigma2=1){mses(n, sigma2)/msesighat(n, sigma2)}
n <- c(2:10,25,100,1000,10000,100000)
sigma2 <- rep(1, length(n))
data.frame(n, sigma2, bs=bs(n, sigma2), bsighat=bsighat(n, sigma2), vars= vars(n, sigma2), varsighat=varsighat(n, sigma2), mses = mses(n, sigma2), msesighat = msesighat(n, sigma2), re = re(n, sigma2))## n sigma2 bs bsighat vars varsighat mses
## 1 2 1 0 0.5000000 2.0000000000 0.5000000000 2.0000000000
## 2 3 1 0 0.3333333 1.0000000000 0.4444444444 1.0000000000
## 3 4 1 0 0.2500000 0.6666666667 0.3750000000 0.6666666667
## 4 5 1 0 0.2000000 0.5000000000 0.3200000000 0.5000000000
## 5 6 1 0 0.1666667 0.4000000000 0.2777777778 0.4000000000
## 6 7 1 0 0.1428571 0.3333333333 0.2448979592 0.3333333333
## 7 8 1 0 0.1250000 0.2857142857 0.2187500000 0.2857142857
## 8 9 1 0 0.1111111 0.2500000000 0.1975308642 0.2500000000
## 9 10 1 0 0.1000000 0.2222222222 0.1800000000 0.2222222222
## 10 25 1 0 0.0400000 0.0833333333 0.0768000000 0.0833333333
## 11 100 1 0 0.0100000 0.0202020202 0.0198000000 0.0202020202
## 12 1000 1 0 0.0010000 0.0020020020 0.0019980000 0.0020020020
## 13 10000 1 0 0.0001000 0.0002000200 0.0001999800 0.0002000200
## 14 100000 1 0 0.0000100 0.0000200002 0.0000199998 0.0000200002
## msesighat re
## 1 0.7500000000 2.666667
## 2 0.5555555556 1.800000
## 3 0.4375000000 1.523810
## 4 0.3600000000 1.388889
## 5 0.3055555556 1.309091
## 6 0.2653061224 1.256410
## 7 0.2343750000 1.219048
## 8 0.2098765432 1.191176
## 9 0.1900000000 1.169591
## 10 0.0784000000 1.062925
## 11 0.0199000000 1.015177
## 12 0.0019990000 1.001502
## 13 0.0001999900 1.000150
## 14 0.0000199999 1.000015 n <- 2:100
sigma2 <- rep(1, length(n))
plot(n, mses(n, sigma2), type="l", lty=1, ylab="mse, re")
points(n, msesighat(n, sigma2), type="l", lty=2 )
points(n, re(n, sigma2), type="l", col="red")
abline(h=1, lty=3, col="green")
n <- 2:100
sigma2 <- rep(10, length(n))
plot(n, re(n, sigma2), type="l", lty=1, ylab="re")
abline(h=1, lty=3, col="green")